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We present a numerical method based on real-space renormalization that outputs the exact ground space of 
"frustration-free" Hamiltonians. The complexity of our method is polynomial in the degeneracy of the ground 
spaces of the Hamiltonians involved in the renormalization steps. We apply the method to obtain the full ground 
spaces of two spin systems. The first system is a spin- 1/2 Heisenberg model with four-spin cyclic-exchange 
interactions defined on a square lattice. In this case, we study finite lattices of up to 160 spins and find a triplet 
ground state that differs from the singlet ground states obtained in CD. Batista and S. Trugman, Phys. Rev. 
Lett. 93, 217202 (2004). We characterize such a triplet state as consisting of a triplon that propagates in a 
background of fluctuating singlet dimers. The second system is a family of spin- 1/2 Heisenberg chains with 
uniaxial exchange anisotropy and next-nearest neighbor interactions. In this case, the method finds a ground- 
space degeneracy that scales quadratically with the system size and outputs the full ground space efficiently. Our 
method can substantially outperform methods based on exact diagonalization and is more efficient than other 
renormalization methods when the ground-space degeneracy is large. 



I. INTRODUCTION 

Renormalization methods are powerful tools for studying 
the long-wavelength properties of physical systems by a sys- 
tematic elimination of high-energy degrees of freedom. The 
first numerical renormalization group (NRG) method was de- 
veloped by Wilson JT] [2| to solve the Kondo problem, an im- 
portant problem in physics that involves the interaction of a 
magnetic impurity with a conduction band (3j. The more re- 
cent density-matrix renormalization group (DMRG) method 
was successfully applied to a large class of one-dimensional 
(D = 1) quantum systems [4| and a few D = 2 sys- 
tems [ 5 - 8 1 . Recent advances in quantum information theory 
also led to renormalization and variational methods, includ- 
ing PEPS mHOl, MERA |JTJ[T2), and tensor renormalization 
lfT3HT6l . The problem with known renormalization methods 
is that they suffer from important limitations when studying 
systems in space dimension D > 2 or with a large number 
of ground states. Our goal is to construct a renormalization 
method that can be applied to such systems when the Hamil- 
tonians under consideration satisfy a "frustration-free" (FF) 
property. 

The term "frustration free" was first coined by the quantum- 
information community [ 1 8 -20 ] to denote a class of Hamil- 
tonians H = X}fe=i n vi, whose ground states are also ground 
states of each local term ir Vk . Vk refers to a finite set of degrees 
of freedom, e.g., a unit or a finite subsystem. While describ- 
ing such Hamiltonians as FF is adequate from a viewpoint that 
we discuss below, the term FF can be confusing if we adopt 
a more traditional convention of identifying frustration with 
competing interactions. For example, the triangular lattice 
Ising model with antiferromagnetic (AFM) exchange (J > 
0) is the paradigmatic example of a frustrated Hamiltonian. 
However, this model is FF according to the previous defini- 
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tion. We let be the three spins = {—1, 1} in the fcth tri- 
angle, j = {1, 2, 3}, and define n Vh = (J/2)[(£j <y)f - 1]. 
The Hamiltonian is FF because any ground state \tp) of H 
satisfies TT Vk \tp) = 0, i.e., \ip) is also a ground state of each 
TT Vk . Nevertheless, does not minimize each of the bond 
Hamiltonians Ja^cr^,, reason why the model is considered to 
be frustrated according to the traditional convention. 

The previous discussion implies that the concept of frustra- 
tion is relative to a particular decomposition of H. The tra- 
ditional interpretation of frustration assumes a decomposition 
of H dictated by the physical nature of the interactions. How- 
ever, while H may be frustrated with respect to one decom- 
position, it may still be FF because the competition between 
interactions on different units disappears when we consider a 
different decomposition (e.g., triangles instead of bonds in the 
Ising example). Remarkably, FF Hamiltonians are ubiquitous 
in condensed matter and quantum information theory. They 
include Ising models, the AKLT model 11211 . parent Hamilto- 
nians of PEPS [9 1, and Hamiltonians that can simulate quan- 
tum circuits [22]. Several other frustrated magnets also corre- 
spond to FF Hamiltonians [23 1. Ground states of FF Hamilto- 
nians contain all the characteristics of highly frustrated phys- 
ical systems: large ground state degeneracy [24 1, coexistence 
of different phases, and exotic orderings. 

In this manuscript, we introduce an exact real-space renor- 
malization method (ERM) that obtains the full ground-space 
of FF Hamiltonians. The output of the ERM is a sequence 
of tensors whose contraction allows us to compute expecta- 
tion values of observables and amplitudes of the ground states 
(Sec. |ll]>. The computational cost of our method (i.e., the cost 
of the tensor contraction) is polynomial in the ground-space 
degeneracy of the FF Hamiltonians involved in the renormal- 
ization steps. If such a degeneracy increases polynomially 
with the system size, the ERM is efficient. Otherwise, for 
exponentially large degeneracies, the ERM is inefficient but 
can substantially outperform other numerical techniques for 
this problem. 

To illustrate the potential of our method, we apply it to two 
FF spin systems that have largely degenerate ground states. 
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The first system is a spin- 1/2 Heisenberg model with four- 
spin cyclic-exchange interactions that is defined on a square 
lattice (Sec. Ill I. The ground state degeneracy is exponen- 
tial in the linear size, L, of the lattice. Such a degeneracy is 

T 2 

much smaller than the Hilbert space dimension 2 , so the 
ERM outperforms exact diagonalization in this case. Besides 
the singlet ground states that were in identified in Ref. 11251 . 
we find a triplet ground state that consists of a triplon that 
propagates in a background of fluctuating singlet dimers. The 
propagation of this triplon leads to an incipient long range 
AFM ordering, which indicates that the triplet ground state de- 
scribes an AFM quantum critical point. Such a triplet ground 
state also exists in rectangular spin lattices of size L x x L y , 
L x < L y , with periodic boundary conditions. In particular, 
Lajko, Sindzingre, and Penc, gave an analytical expression of 
this state for the case of three-leg tubes (L x = 3) [26|. 

The second of system consists of the family of spin- 
1/2 Heisenberg chains with uniaxial exchange anisotropy 
and nearest and next-nearest-neighbor exchange interactions 



(Sec. IV I. The Hamiltonians in this family are also FF. An an- 
alytical and closed form representation for the ground states 
of this family, in terms of anyonic operators, was given in 
Ref. l23ll . However, such a representation may not be useful 
for computing some expectation values of spin-spin correla- 
tions efficiently. Those expectation values can be efficiently 
computed with the ERM. 

We note that quantum Monte Carlo methods, that are not 
based on renormalization, cannot be applied to most FF 
Hamiltonians because of the infamous sign problem. 
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FIG. 1. Spin system. Black dots denote spins located at the vertices 
V of the lattice. Each term 7r„ fe in the Hamiltonian acts nontrivially 
in Vk only. The sets wi and zi are defined in the text. 



The Hamiltonian is represented as 

p 

H = ^ TT Vk 



(1) 



fc=l 



where each n Vk is a Hermitian operator acting nontrivially on 
spins in v k only. We assume ir Vk > 0. If any ground state \ip) 
of H satisfies 



H\^)=ir Vl \^) 







(2) 



II. THE RENORMALIZATION METHOD 

We start by providing a brief description of the ERM (tech- 
nical details are provided in Appendix [A|. For simplicity, we 
consider a Hamiltonian H acting on a system of L spins lo- 
cated on the vertices V of a lattice [30]. The Hilbert space of 
the system is T~L\> and its dimension is dy. The magnitude of 
the spins and space dimensionality of the lattice are arbitrary; 
we refer to the spin system of Fig.[TJfor illustration purposes. 
We let Vk be a subset of spins of V, Hk the associated Hilbert 

k = 1,2, ... ,p and the 
. . , p we also define the 



space, and d Vk its dimension. Here, 
subsets Vk are known. For I = 1,2,. 
sets of spins wi = Vi\zi-±, with zi 



U 



fc=i 



Vk 



U k=1 Wk 



and zq = {0}. A\ B is the relative complement of A in B so 
that wi are those spins that belong to vi but do not belong to 
any other Vk with 1 < k < I — 1. We assume it); ^ {0} and 
note that V\ = Wi = z\. 

For a set x of spins in V, we use i x for the spin variable 
in some standard basis ll32l . Also, d x denotes the dimen- 
sion of H x , the Hilbert space associated with x. It follows 
that {|i x )}i<^<d m is an orthonormal basis for the spins in x; 
hereafter referred to as the computational basis and \i x ) is a 
basis state. In some cases, we do not distinguish between basis 
states or (column) vectors: \i x ) can also denote a vector with 
component equal to 1 in position i x and zeroes elsewhere. The 
number of components of \i x ) is the number of values that i x 
can take (< d x ), which is given in each case. 



then H is said to be FF. The standard definition of FF also 
assumes that each iv Vk is local. Here, we can relax such an as- 
sumption without incurring in large computational overheads 
as long as ir Vk is a bounded sum of product operators (see 
Appendix [A|>. When H is FF, our renormalization method is 
exact and outputs the ground states of H as 



v \ _ 



• T tV v ) 



(3) 



T{ , Tj, . . . , are isometries, 1 < i v < g, and g is the 
ground-space dimension. The symbol • refers to a tensor con- 
traction that involves a sum over repeated indices. Such a con- 
traction regards a tree-tensor network (see Fig.|2]i; Tree-tensor 
networks are ubiquitous in renormalization methods Il33ll34l . 
Whether H is FF or not is also an output of the ERM. 

The ERM performs p steps. Each step can be defined re- 
cursively as follows. For 1 < I < p, we let g\ > be the 

ground-space degeneracy of Hi = J\=i n vh (<?o = !)■ The 
ground states of Hi are 1 < i Zl < gi- H p = H. In 

the first step, the ERM diagonalizes 71, the d Wl x d Wl matrix 
representation of ir Vl in the computational basis. It obtains 
gi and {l^i )}i<^i<si' an d continues only if 171 > 0. In 
the Z-th step, / > 2, the ERM computes 7; . This is a hi x hi 
matrix representation of ir Vl in the basis , i Wl }}, with 

1 < ^'-1 < 9l _ u 1 < i w ' < d Wl , and h t = gi-id wy 
Then ERM applies exact diagonalization to ji, obtains a basis 
{\iplz t )} f° r the zero eigenvalue, and computes the multiplic- 
ity gi of the new ground space. The ERM continues only if 
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FIG. 2. Representation of the tree-tensor network contraction for a 
system with p spins (blue circles). In this example, the sets vi refer 
to pairs of nearest-neighbor spins and each wi is a single spin. The 
tensors T; live in the vertices of the tree and the contraction indices 
are in the edges. An arrow means a sum over the corresponding 
index. The contraction T p • . . . • Ti maps ground states of H into 
states |i z p).T 1 t »---»T p t isthe inverse transformation that gives all 
the ground states of H from \i Zp ). 

gi > 0. The isometries in Eq. ^ are 

91 

t}= £ lvi-«)<i*M- w 

In Appendix |A] we show that, if Nt is the number of ele- 
mentary operations to output the isometries T± , . . . ,T£, then 

p 

N T cxY,EV(h k ) + (p-k)(h k g k ) 2 . (5) 

EV(d) is the cost of the exact diagonalization of a dx d matrix 
and go = 1. Also, if Mt is the memory cost associated with 
the number of variables kept during the implementation of the 
ERM, 

p 

M T oc ^2 h k9k ■ (6) 

k=l 

Thus, the efficiency of the ERM strongly depends on g k and 
the method becomes efficient when g k £ 0[poly(p)]. The 
cost of evaluating expectation values of observables in any 
ground state of H is also important and can be easily derived 
from the analysis given in Appendix [A] 

A related method for solving some spin- 1/2 systems, which 
is based on the techniques developed in Ref. |27l . can be 
found in Ref. l28l . Also, an exact renormalization method 
for quantum spin chains was proposed in Ref. [29]. Our 
main contribution with respect to that of Refs. Il28l |29l is that 
we consider Hamiltonians with arbitrary interactions, in any 
space dimension, and provide a real-space renormalization al- 
gorithm that is exact if the ground space satisfies some proper- 
ties that can be verified by the ERM. In addition, the way that 
the ERM contracts tensors is different from the usual contrac- 
tion of a binary-tree like tensor network. The main reason be- 
hind this difference is the minimization of Mx when gi grows 
monotonically with I. 



III. APPLICATION TO A D = 2 MAGNET 

A. Model Hamiltonian 

We apply the ERM to a spin- 1/2 model that satisfies the 
frustration-free property as defined in Sec. [II] The FF Hamil- 
tonian is defined on a square lattice and reads [25]: 

Each operator V a projects the total spin state of a square 
plaquette a onto the subspace with spin 2, as illustrated in 
Fig(3ja). This model has a number of exact valence-bond or- 
dered ground states that increases exponentially in the linear 
dimension of the square lattice 11251 . The model of Eq. |7]) 
corresponds to a particular regime of parameters of a more 
general model with frustration and ring exchange: 

H' = J\ s r • s r / + J 2 ^ s r ■ s r / 

(r,r'> «r,r'» 

+ K (P"Pki + PfuPa +P? k P?i)- (8) 

Here, (r, r') and ((r, r')) denote nearest neighbors and next 
nearest neighbors, respectively, and i, j, k, and I label the 
four spins of a plaquette in cyclic order. s r = (s^ , s;f, s^) 
is the spin operator of the spin at the rth position. Equation 
|7]l reduces to Eq. up to an irrelevant constant, if J\ — 
1, J2 = 1/2, and K = 1/8. This case corresponds to the 
maximally frustrated point as a function of J2 / J\ ■ 

B. Algorithm 

The implementation of the ERM is related to Wilson's NRG 
[I] |3 and to the warmup stage of the conventional DMRG 
Thus, it is simple to adapt an existing NRG or DMRG 
code for this case by making minor modifications. At each 
renormalization step, the ERM grows the system by adding 
one plaquette. Then, the ERM diagonalizes the renormalized 
Hamiltonian and only keeps the ground states (if the lowest 
eigenvalue is zero). In a DMRG "language", this corresponds 
to working with two blocks instead of four. 

The order in which the plaquettes are added can be arbi- 
trary. However, the ground space degeneracy, <?;, depends 
dramatically on the path followed to grow the lattice. The 
"snake" path shown in Fig. |3jb) is the approach that turned 
out to be more efficient. The number of spins increases by 
one or by two at each step [see for instance steps 5 and 6 in 
Fig. (3jb)]. It is important to introduce Hamiltonian terms cor- 
responding to a plaquette at a time (the projectors V a ) and 
to avoid including terms belonging to neighboring plaquettes. 
For instance, a nearest neighbor exchange term in Eq. ([8| is 
shared by two neighboring plaquettes. Such a term should be 
split accordingly to assure that only one projector is added per 
step. 

We applied the ERM to rectangular lattices of sizes L x x 
L y , L x < Ly, and periodic boundary conditions. The ground 
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FIG. 3. (Color online) (a) Square lattice and the projectors V a acting on the square plaquettes (dashed lines). Each V a projects the spin states 
of the plaquette onto the subspace with spin 2. (b) "Snake" path followed by the ERM, similar to the DMRG warmup scheme. At each step 
one adds the spins necessary to complete a projector on a new plaquette. This growing process requires adding one or two spins (thick bold 
lines) per step. 



space degeneracy g is proportional to 2 Lm in this case, i.e. ex- 
ponential in the linear dimension. If we follow a path like the 
"snake" in Fig. |3jb), gi increases as we increase L x but it does 
not change substantially when we increase L y . This property 
allows us to study remarkably large system sizes by keeping 
L x constant and by increasing the number of plaquettes in 
the y direction to relatively large values of L y . Every time 
we close a boundary along the y direction, the degeneracy gi 
drops substantially. 

Unlike DMRG, the ERM needs to fully diagonalize the 
renormalized Hamiltonian at each step. This implies diago- 
nalizing a hi x hi matrix in the Zth step, with hi = 2gi_i 
or hi = ^gi-i, depending on the number of spins added at 
that step (one or two, respectively). The maximum value of 
hi depends on the linear dimensions of the system, and could 
reach several thousands for lattices of hundreds of spins. To 
be able to study relatively large systems, we can use symme- 
tries -U(l)/abelian quantum numbers in our case- to store the 
Hamiltonians and other operators in block form. However, to 
ensure that we keep all the ground states, we do not restrict 
the values of these quantum numbers. 

In Fig. |4]we plot the ground-space degeneracy, gi, and or- 
der of the renormalized Hamiltonian, hi, for a 6 x 6 lattice 
(L x = L y = 6). We also plot the order of the largest block 
of the renormalized Hamiltonian when symmetries are con- 
sidered, ft," lax , which determines the dominating cost of the 
method. The oscillation in gi has a periodicity corresponding 
to the linear dimension of the lattice, showing the reduction in 
the ground space degeneracy every time a lattice boundary is 
closed. 

Because we are interested in the computation of correlation 
functions of arbitrary ground states, at each step we need to 
store all the matrices of the operators involved in such cor- 
relations. We note that correlators of the form (AB) cannot 
be computed by storing the operators A and B independently. 
In general, it is necessary to store the matrices for the prod- 
uct AB because the product of two projected operators is not 
equal to the projection of the product. 

The ground states of H can be obtained for 8 x 8 and larger 




step I 

FIG. 4. Ground space degeneracy (<?;), order of the renormalized 
Hamiltonian (hi), and order of the largest block in the renormalized 
Hamiltonian (ft™ ax ) per step of the ERM, when symmetries are con- 
sidered. 

lattices by implementing the ERM on more powerful exist- 
ing computers. The requirement of storing all the operators 
associated with the correlations described in the next section 
and those needed for computing the Hamiltonian terms, is the 
main limiting factor for increasing the lattice size. 

C. Results 

In addition to the large set of exact 5 = ground states ex- 
hibiting valence bond ordering [25], the results output by the 
ERM show a triplet ground state, with S = 1, S z = ±1, 0. (S 
is the total spin of the lattice and S z the zth component of the 
total spin.) Such a state was also identified by exact diagonal - 
ization of small square clusters and by an analytical solution 
of the model on 3 x L y tubes |26l[35). That the S = 1 state 
exists in larger systems was unexpected and illustrates the im- 
portance of methods that obtain the full ground space. The 
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existence of a triplet ground state is compatible with a critical 
scenario in which H has a gapless spectrum of S = 1 spin 
excitations. In this case, an external magnetic field B > 
would induce AFM ordering of the spin components that are 
orthogonal to the field's direction. This scenario is confirmed 
by the spin-spin correlators obtained with our ERM for the 
S = S z = 1 ground state on 4 x L y and 6 x L y finite lat- 
tices. Figure [5] shows the two-point correlators (sqS*) and 
(sqS~) as a function of the distance along the y-direction. 
= (0, 0) denotes a reference spin (the origin), r = (0, y), 
and sf = ± is^.. While there is a clear long-range AFM 
tail for (sqS~), the (s s r) correlator decays exponentially 
in r. It is important to note that the magnitude of the local 
staggered a;y-magnetization is of order 1 / yjL x L y , which is 
the expected behavior for the condensation of a single triplon. 
This property is revealed by the scaling behavior of the struc- 
ture factor, 



S ± (k) 



1 

L x L n , 



E 



eJ k.(r-r') (s + s - ) 



(9) 



evaluated at the AFM wave vector k = (tt, tt) (inset of Fig.[5]l. 
S ± (tt, tt) tends to a value of order one for L x , L y — »• oo, in- 
dicating that the the order parameter L l L ^ r e lk r s^, v — 

{x, y}, is proportional to \j JL x L y . 
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It can either select one particular valence bond ordering via 
an order by disorder mechanism, or simply destroy any long 
range bond ordering. In the former scenario, the application 
of a small magnetic field, B, that couples to the spins via the 
Zeeman term — B^ r S z , should stabilize a particular bond 
ordering out of the exponentially large number of degenerate 
bond ordered ground states that exist at B = 0. The selection 
mechanism would be provided by the kinetic energy of the 
field induced triplons, which should be minimized for a par- 
ticular bond ordered background. In this scenario, the selected 
bond ordering coexists with the AFM a;y-ordering induced by 
the condensation of triplons at the single particle state with 
momentum k = (tt,tt). Figure [6] illustrates this situation for 
the bond ordering that has the highest susceptibility, according 
to the results that we discuss below. In the second scenario, 
the bond fluctuations induced by the triplon propagation are 
strong enough to produce a valence bond liquid with short- 
range bond-bond correlations. The only order parameter that 
survives is the AFM ordering which arises from the triplon 
condensation. 

To further explore both scenarios, it is necessary to compute 
bond-bond correlation functions for the S — S z — 1 ground 
state. We introduce the bond structure factor, Sb (k, v), for the 
local bond operator B rl/ = s r s r+e ^ . k = (k x , k y ) is the wave 
vector, v — {x, y}, and e„ is the relative vector connecting 
nearest-neighbor spins along the ^-direction. Then, 

S*(k,i/) = — L- Ve lk '( r - r ')(i? r ^ r v) ■ (10) 

L x -L y 

y r,r' 

The staggered bond ordering (SBO) shown in Fig. [6] should 
produce a sharp maximum in 5s(k, v) at k = (n, tt). In 
contrast, the bond structure factors obtained for 4 x 4, 6 x 6 
lattices have very broad maxima at k = (tt, tt), indicating 
that the second scenario with short ranged bond correlations 
is more appropriate for the (S — l,S z — 1) ground state 
(see Fig. [7] a and b). The same scenario hods for the 4 x 40 
lattice (Fig.|7]c). In this case, Sb (k, y) has an approximately 
degenerate line of maxima for k x = tt/2, indicating that bond 
correlations between adjacent vertical lines are very weak. 



IV. APPLICATION TO A D = 1 ANISOTROPIC 
HEISENBERG MODEL WITH NNN INTERACTIONS 

The one-dimensional family of spin- 1/2 Hamiltonians in- 
troduced in Refs. [361 and [23] is 



FIG. 5. Two-point correlators {sgsl) and (sq s^) as a function of 
distance along the y-direction. Here, is a spin of reference and 
r = (0, y). The figure shows the correlations for a) 4 x L y and b) 
6 x L y lattices. Inset of a): scaling of the structure factor, S ± (k), 
evaluated at the AFM wave vector k = (tt, tt). 

This incipient AFM ordering occurs at the same point 
where an exponentially large number of different valence or- 
derings become degenerate (5 = ground state sector). How- 
ever, the presence of a triplon that propagates across the lattice 
could lead to two different scenarios for the bond correlations. 



J>=1,2 



A, 



a x x i y y 



The Hamiltonian coefficients are parametrized by the single 
variable Q: A„ = cos(isQ) and J\ — —4J2 cosQ, with 
< Q < tt. Hq satisfies the FF property for all Q if we 
choose appropriate boundary conditions. In particular, for 
Q = tt/2, the model reduces to two decoupled ferromagnetic 
Heisenberg chains whose exact solutions are known (37). In 
this case, the ground-space dimension is g = (L/2 + l) 2 
or g = (L + 3) (L + l)/4 for even or odd L, respectively. 
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FIG. 6. Illustration of one of the possible scenarios in presence of 
a finite magnetic field B. The ovals indicate that the corresponding 
bonds have a predominant singlet character. The staggered bond or- 
dering (broken Z4 symmetry) shown in the figure should produce a 
sharp maximum in the bond structure factor Sb (k, v) at k = (tt, 7t). 
The arrows indicate the AFM ordering in the plane perpendicular to 
the applied magnetic field induced by the condensation of triplons. 
The staggered magnetization can point along any direction obtained 
by a global spin rotation along the field direction (broken U(l) sym- 
metry). 



L is the number of spins in the chain. The analytical solu- 
tions for the ground states presented in Ref. E3ll show that the 
ground space for Q = ir/2 is continuously connected with the 
ground space for arbitrary Q. While, in principle, the ground- 
space degeneracy for arbitrary Q could be larger than that for 
Q = tt/2, the ERM shows that such a degeneracy remains 
constant. Therefore, the ERM allows for an efficient compu- 
tation of spin-spin correlations in any ground state of Hq. 



V. CONCLUSIONS 

We introduced an exact renormalization method that ob- 
tains the full ground-space of Hamiltonians satisfying a 
frustration-free property. The method outputs a sequence of 
tensors whose contraction allows for the computation of cor- 
relation functions in any ground state. Such correlations can 
be used to characterize zero-temperature states of matter. The 
cost of the method depends on the ground state degeneracy for 
each renormalization step. The method computes the degen- 
eracy and verifies the frustration-free property. 

We applied the method to two spin systems. First, we con- 
sidered a FF Hamiltonian in D — 2 whose ground-space di- 
mension increases exponentially in the linear size of a square 
lattice. We applied the ERM successfully and characterized 
the physical properties of the only ground state with total spin 
5 = 1 and projection S z = 1. Such a ground state differs 
qualitatively from the (bond ordered) singlet ground states 
identified in Ref. Il25l . In particular, it describes a quantum 
critical point associated with the onset of AFM ordering. Ac- 
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FIG. 7. Bond structure factors Ssf^k, x) as a function of k = 



, computed with the ERM for different lattice sizes, a) L x 



4. b) L x = L y = 6. c) L x = 4, L y 



40. 
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cording to our results, the application of an arbitrary small 
magnetic field B > induces AFM order of the spin com- 
ponents that are orthogonal to field's direction. Our results 
also indicate that the finite magnetic field destroys any of the 
long-range bond orderings that compete at B = l25ll . 

We also applied the ERM to the one-dimensional family of 
FF spin- 1/2 Hamiltonians described in Eq. ( [TT| . We verified 
that the ground space dimension is g = (L/2 + l) 2 (g = 
(L+3) (i+l)/4) for even (odd) L, as conjectured in Ref. Il23l . 
The ERM is efficient in this case. 

Both applications illustrate the power of the ERM for ob- 
taining exact ground state properties of frustration free Hamil- 
tonians. Nevertheless, our renormalization method can also be 
used in more general contexts. For example, since stochastic 
matrices can be sometimes related to FF Hamiltonians ll38l . 
the ERM can be used to compute properties of (classical) sys- 
tems in and out of equilibrium. Similarly, the ERM can be 
used to solve combinatorial optimization problems [39|. 
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Appendix A: Exactness and cost of the ERM 

We provide more details about the ERM, showing the ex- 
actness of the method and estimating the computational re- 
quirements. With no loss of generality ir Vk is 

^=e°m (i) ®-® o ;- w . < ai > 

a=l 

where ® is the tensor product. The sets low^) satisfy W^ij) H 
v k 7^ {0}- Any state \i Vk ) corresponds to a particular state 
\i Wkl - r > , . . . , i Wk (<~)) for the sets w k (j) in Vk- We assume then a 
specification of H via access to that, on input (v^, i, i'), 

it outputs all the matrix elements 

(i' Wk ^\0^\t Wk ^) , (A2) 

for 1 < j < r and 1 < aAk- ITh also outputs all the sets 
Wyj) involved in each term of the decomposition of ir Vk . In 
other words, Hh gives all the information about the action of 
each term of ir Vk in the states \i Vk ). Such a Hamiltonian spec- 
ification is common in applications, e.g. a spin- 1/2 system 
specified in terms of Pauli operators. The actual computa- 
tional cost of the ERM should consider the number of times 
that n# is used, which is typically linear in p. 

Our renormalization method takes the ordered set 
{i>i, . . . Vp] as input and uses Hh- It outputs a bit b denot- 
ing whether H is frustration free (b = 0) or not (b = 1) ll40l . 
If b = 0, the ERM also outputs the ground states specified by 



a sequence of tensors 7] or Tj - see Eq. Q. Because Eq. ([3]) 
involves a summation over repeated indexes, we use Einstein 
notation in the following. For / = 1, . . . ,p, T\ is determined 
by its (complex) entries: 

f if/ = l 

Ti:={ , (A3) 

if/>i 

with 1 < i Wl < d Wl , 1 < i Zl < gi, and gi determined by 
the ERM (see below). For each i Zl < gi, the vectors ) = 
Ti»...»T}\i Zl ) arein% 2l and thus haved 2i components in the 
computational basis. We recall that zi = U^. =1 Wfc. Then, each 
such component is determined by i Wl , . . . , i Wl corresponding 
to the sets w\, . . . , wi, respectively. That is, {\i Wl , . . . ,i Wl )}, 
with 1 < i Wk < d Wk , also defines a computational basis for 
H Zl ■ The components of \ip\ H ) in such a basis are 

(u^Qiu^ 2 ...(u^- 1 ^ 1 , (A4) 

with 

(ui)Q =((*i)A)* (A5) 

(u^tT 1 '^" ={{t k %u^ k y , 2<k< P . 

In particular, when I — p, we have z p = V and Eq. \AA\ de- 
fines the contraction in Eq. (j3J. This contraction is associated 
with a tree-like tensor network - see Fig[2]for an example. 

For I = 1, we denote by 71 the d Vl -dimensional matrix rep- 
resentation of 7r„j in H Vl — H Wl , in the computational basis. 
The ERM constructs 71 using n# once. Then, the ERM per- 
forms exact diagonalization to obtain g\ and an orthonormal 
vector basis {\4>}z 1 )}i<i^i< gi for the zero-eigenvalue eigen- 
vectors of 71. T\ is the tensor that maps such eigenvectors to 
vectors or states in the computational basis: 

91 

Ti= N Zl )(0kl • (A6) 
1-1=1 

The entries of T\ are 

(ti)& =(4>U\i Wl ) » (A7) 

with 1 < i Wl < d Wl and 1 < z x < g x . 

In the Z-th step, I > 2, the ERM uses to construct the 
hi x hi matrix 7/, hi = gi-i-d Wl , with entries 

Jl^i^L,^ 1 ]^ 1 ^,^ 1 ) ■ (A8) 

Here, = 1, . . . , g l _ 1 and i' Wl , i w ' = l,...,d wr 

The ERM performs exact diagonalization of 7; and continues 
only if the lowest eigenvalue is (b = 0). It computes the mul- 
tiplicity of the zero eigenvalue, assigns it to gi, and computes 
a complete orthonormal basis of hi -dimensional eigenvectors 
■ ■ ■ , \4> l gi )} f° r me zero eigenvalue. The ERM assigns 
the tensor T ; to the transformation that maps such eigenvec- 
tors to vectors in the computational basis of % ZI : 

91 

T l =Y i \i*)<A-A- (A9) 
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The entries of 7] are 



(A10) 



To show that the ERM is exact, we first remark that the 
eigenvectors \4>\n) represent the states \ip^i), as determined 
by Eq. \AA\ . Then, we will show that the states 



\4 



Tl»...»T?\i l 



l \"zl ' 



(All) 



given according to Eq. ( |A4| i, are ground states of Hi = 
Sfc=i n vk f° r a ll I € {1,2,... when 7J/ are frustration 
free. 

The proof is inductive. For / = 1, \<f>\z 1 ) = is a 

ground state of jx or 7r Vl = i?i by definition. That is, ), 
as determined from Eq. ( |A4| i, is a ground state of Hi for each 
z Zl = 1, . . . , g\. We assume now that 



{|^ 1 > = T 1 t .....T/_ 1 i4- 1 » 1 < i , ! _ 1 



<gi-i 



is an orthogonal basis for the ground subspace of Then, 
if |0) is a ground state of Hi, 



gi-i a ^i 

E E 



(A12) 



withc^i-i ^w, complex amplitudes. Otherwise, |0) would not 
belong to the intersection between the ground subspaces of 
Hi-x and Hi, a requirement for frustration-free Hamiltonians. 
It suffices to obtain 7;, a projection of tt Vi into the subspace 
spanned by , i*")} i<i*'-i< S! _i ) i<i'»!<(i, ■ Without 

loss of generality, we can choose a value of i zi such that 



Thus, 



\4>) = 14 



(A14) 



where ' are the entries of Tj . That is, 

\<t>) = {ui%T 1 ^ l We l -^ Wl ) 



(A15) 



(ui)^ 1 (u 2 ) • 



The contraction in Eq. ( |A15[ > coincides with that of Eq. ( |A4| i. 
It follows that 



lW-,>=7? 



• T, T li 2 



(A16) 



is a ground state of Hi for all 1 < i Zl < gi- In particular, 
\if>f'p), with 1 < i Zp < g p = g, are all the ground states of H. 



the memory requirements, i.e., the number of coefficients that 
need to be kept in memory during the implementation of he 
ERM. 

We first write 



N 7 



E^ 



V 

1=1 



(A17) 



To obtain T\, the ERM performs exact diagonalization of 
71. Because 71 is of dimension d Wl , the cost of obtain- 
ing its eigenvectors and eigenvalues is TYj. < EV(d Wl ) g 
0(poly(d Wl )). E'y(d) is the cost of exact diagonalization of 
a d x d matrix, which is almost quadratic in d in actual im- 
plementations. Only the g\ eigenvectors with zero eigenvalue 
need to be kept in memory for the following step and thus 
Mj. ex g\d Wl . 

The cost of obtaining all zero-eigenvalue eigenvectors of 7; 
is bounded by EV(hi). To obtain Nj, we need to add the cost 
of computing 7; . Each matrix element of 7/ in Eq. dA8) is 



(i/ !S '-i\(i; w >\T l _ 1 .....T 1 .ir Vl . 
•T± •...•T/_ 1 |i*<- 1 )|^<) . 



(A18) 



We consider the decomposition in Eq. jAl\ and we are inter- 
ested in obtaining the cost of computing a particular term 



(i' z ^\{i ,Wl \T l _ 1 .....T 1 .(o^®.. 

...®or,r).Tt.....Tt_ 1 |^-)r) 



(A19) 



We can interleave trivial operators t w = 53i«=i \i w )(i w \ m 
Eq. ( |A19| l for those w 7^ Wk(j) without affecting the output of 



(A13) the contraction. That is, we extend the definition of O 



WktJ 



that 0^ <J) = liufcy) f° r r < 3 ' < Then we write 
for all 1 < j < I. Equation ( |A19| l is 



(top.-.ftj-oj:::: 



(A20) 



(A21) 



{oita («i)W-(«j-o;.«:x 



z i-2_,;™i-i 



where we used Eq. ( |A4| >. As before, Eq. ( |A21[ ) refers to a 
contraction of a tree-like tensor network; see Fig. [8] for an 
example. 

The cost of evaluating Eq. ( |A21[ > depends on the support 
of vi, that is, the number and position of spins that belong to 
vi. This is so because Eq. ( |A2 1 [ > can be sometimes simplified 
considering that 



1. Computational requirements 



We let Nt be the total cost, i.e., the total number of elemen- 
tary operations to obtain all entries of T\, . . . ,T p . For sim- 
plicity, we do not consider in the cost the number of queries 
to Uh, which is typically linear in p. We also let Mr be 



and then 



9k 



T k Tt= £ |i**><»* 



(A22) 



i« h 1,1 k , (A23) 
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FIG. 8. Representation of the network contraction for each ma- 
trix element of 7; in Eq. «[A2TJ for the same system of Fig. [2] 
Blue circles are spins. n Vl is a sum of At terms of the form 
O™^ Cg> . . . <S> 0™ l a , and some O k 3 a may act trivially in Wj\ i.e., 

Ofea = S •«"/_! |* u '-'}(*"^|- Arrows denote a sum of the index in 
the corresponding edge. Open circles denote a fixed index, refer- 
ring to a particular matrix element of 71, the projection of n vi in the 
ground subspace of H1-1 = X/fc=i 1 — i Wl ,i' Wl < du^ and 

1 < i'* 1 - 1 < 

which are useful if 0^ <j) = ^-w kuy Nevertheless, to analyze 
the cost of computing Eq. ( |A21| i we consider the worst case 
scenario in which O k ^ <j) 7^ %n k{j) f° r a ll 1 < 3 < I- We can 
rearrange the sum and compute Eq. ( |A21| > in I sequential steps 
as follows. First, we compute the 

for all 1 < i Zl ,i' Zl < g x . Because 1 < i w \i' Wl < d Wl , this 
step has cost oc (d Wl gi) 2 . We keep the computed values in 
memory. Next we compute the 

(tiJP, (ta)^,^ 

for all 1 < i Z2 ,i' Z2 < g 2 . This step has an additional cost 
oc (gid W2 g2) 2 , where the first gi comes from the sum in i Zl 
and i' Zl . We keep implementing the procedure sequentially 
until we compute 

(^(oi^MP... (A24) 

■ ■ ■ [tl-l) i »t-3 ii H'l-l \V'l-l) i z l _ 1 , 

which has an additional cost oc (/i;_ig;_i) 2 with respect to 
the previous computations. That is, the sequential method has 

an overall cost ex Y^k—ii^k-gk) 2 , w i m 9o = !■ The sequen- 
tial method can be understood from the example in Fig. [8] in 
which the sequential steps regard the contraction of tensors 
from left to right. 

The last step is to compute 

(oi)f-i (A25) 

for all 1 < i Wl , i' Wl < d W[ . This step is implemented us- 
ing Hh and has no cost under our assumption. Then, the 



computation of Eq. < [A2j) for all 1 < i w ' , i' w ' < d Wl and 

1 < i 21 - 1 , i' 21 - 1 < gi_i can be implemented with 

K ^2(9k~id Wk .gk) 2 (A26) 
fc=i 

elementary operations. To compute ji, we need to add a mul- 
tiplicative factor Ai that regards the number of terms in the 
decomposition of 7T„, . Then, 

1-1 

N l T = EV( gi ^d Wl ) + cA t J2( h k-9k) 2 , (A27) 

k=l 

with c > 1 a constant. A memory of M| oc higi is needed for 
the ground states of 7;. 

Typically, Ai is bounded by some constant A. In this case, 

p 

N T oc EV{h k ) + (p- l)(d Wl g l ) 2 + 

k=l 

+ (p - 2)(/ l252 ) 2 + . . . + ( Vi3p-i) 2 ■ (A28) 
In addition, 

p 

M T oc ^ /i fc5fc . (A29) 
fc=i 

If 5 fe € C[poly(p)] and G 0[poly(»], then JVy G 
0[poly(p)] and the ERM is efficient. 

2. Optimal cost 

The total cost of the ERM in Eq. ( |A28| l depends on d Wl and 
gi. In many applications, d Wl is constant and Nt and Mr 
are functions of g x , . . . ,g p . The cost can then be minimized 
by considering all possible orderings of w\, . . . , w p such that 
Nt and/or Mt are minimum. This procedure rules out some 
possible orderings that yield exponential complexity in sys- 
tems in which the ERM could be implemented efficiently. For 
example, consider a square spin lattice and assume that the 
terms n Vk in the Hamiltonians involve two nearest-neighbor 
spins in either direction. The ERM could be implemented to 
obtain the ground states of each chain along a particular direc- 
tion and then add the Hamiltonian terms in the other direction. 
However, this construction results in exponential complexity 
if the ground space of each chain is degenerate because the 
number of ground states that have to be kept in memory is ex- 
ponentially large in the length of the chains. A more efficient 
choice considers "growing" the system using thr "snake" path 
depicted in Fig. [3](b). 

In a different example we consider a binary tree of depth 
q and assume p = 2 q . Each node in the basis of the tree 
corresponds to a single spin in a lattice. A standard real- 
space renormalization method for such a binary tree will have 
/ = l,..,,q steps, each involving a diagonalization of 2 q ~ l 
matrices of dimension g 2 i each. The memory requirement 
for such method is dominated by the last step, which requires 
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dealing with a subspace of dimension g p / 2 x g p /2, spanned 
by all the ground states obtained in the previous step. In ad- 
dition, each such ground state has (gp/4) 2 components in a 
computational basis. If g^ oc k 13 , the memory requirement 



to implement the last step is M' T £ 0[p 413 ]. Nevertheless, our 
ERM implies M T e 0{p 215 ) in this case. Clearly, M' T > M T 
when /? > 0, p 1. The standard renormalization method 
may outperform the ERM only if ft = 0. 
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